adopath ++ D:\Data\WorkData\702092\ado
* need to indicate adopath and have latest David Roodman's --cmp.ado-- command installed

clear all
set mem 5000m
set more off			
capture log close
log using het_IGE.log,replace 

log off

forvalues i=0(1)2 {
	use earnings_long,clear
	quietly summ cohort if bob==`i'
	global min_`i'=_result(5)
	global max_`i'=_result(6)
	forvalues y=${min_`i'}(3)${max_`i'} {
		use earnings_long,clear
		keep if bob == `i'
		ta year, ge(years)
		di `y'
		ge agesq=ageC^2
		regress  logearn years* ageC agesq if cohort==`y',noc
		predict res if cohort==`y',resid
		keep pnr pnrf logearn res year yob bob  age* cohort 
		keep if cohort==`y'
		compress
		save res_`i'_`y',replace
	}
}


forvalues i=0(1)2 {
	use earnings_long,clear
	quietly summ cohort if bob==`i'
	global min_`i'=_result(5)
	global max_`i'=_result(6)
	use res_`i'_${min_`i'}
	local j = ${min_`i'} + 3
	forvalues y = `j'(3)${max_`i'} {			/*puts residualised earnings back together by bob*/
		append using res_`i'_`y'
		erase res_`i'_`y'.dta
	}
	
	bysort pnr: ge mres=sum(res)/_N				/*long term earnings as the average earnings in that range*/
	bysort pnr: keep if _n==_N
	
	reg mres
	predict rmres, resid						/*take deviations from generational means*/
	keep pnrf rmres bob
	save res_`i',replace
	erase  res_`i'_${min_`i'}.dta
}

use res_0
rename rmres rmresf
drop bob
sort pnrf
save,replace

use res_1
append using res_2
sort pnrf
merge pnrf using res_0

keep if _merge == 3



reg rmres rmresf if bob==1						/*IG regressions, by siblings and joining siblings with family-level clusters in s.e.*/
reg rmres rmresf if bob==2
reg rmres rmresf , r cl(pnrf)

log on
nlcom IGE:_b[rmresf]	,post							/* (1) the average IGE in the population*/
log off

mat beta= e(b)
mat vbeta= e(V)

destring pnrf, gen(rino) 

*reml regressions including parental income in the fixed and/or random effect
mixed rmres || rino: , reml			

log on
/*computes the sib corr from the REML*/
nlcom rho: exp(_b[lns1_1_1:_cons]  )^2 / (exp(_b[lns1_1_1:_cons]  )^2 + exp(_b[lnsig_e:_cons]  )^2) , post
log off
 
mat beta= beta,e(b)
mat zero = (0)
mat vbeta= (vbeta,zero)\(zero,e(V))
mat colnames vbeta = IGE rho
mat rownames vbeta = IGE rho


*tell stata which is the right Varcov
cap program drop foo2
program define foo2, eclass
	tempname b V
	matrix `b' = e(b)
	matrix `V' = e(V)
	matrix define `b' = beta
	matrix define `V' = vbeta
	ereturn post `b' `V' 
end

foo2

log on
nlcom  solon_share: _b[IGE]^2/_b[rho]
log off
 
summ rmresf if bob==1,detail
 
sca varX=r(Var)				/*var(y_f)*/ 
 
summ rmres ,detail
 
sca varY=r(Var)				/*var(y)*/ 

log on
*Solon/Bjorklund decomposition without stationarity assumption
nlcom (IGC:			_b[IGE] * sqrt(varX/varY) ) (rho:_b[rho])	,post
nlcom Anders_share: _b[IGC]^2 / _b[rho]							
log off


 
mixed rmres rmresf || rino: , reml  /*this regression applies the sequential conditioning approach,*/ 
									
log on
/*looks at the sibling correlation after including y_f on the rhs of the mean (fixed) effects*/
nlcom rho_bash: exp(_b[lns1_1_1:_cons]  )^2 / (exp(_b[lns1_1_1:_cons]  )^2 + exp(_b[lnsig_e:_cons]  )^2), post
log off
mat beta= beta,e(b)
mat zero = (0,0)
mat vbeta= (vbeta,zero')\(zero,e(V))
mat colnames vbeta = IGE rho rho_bash
mat rownames vbeta = IGE rho rho_bash

foo2
log on 
nlcom bash_share: 1 - _b[rho_bash]/_b[rho]
log off



mixed rmres rmresf ,noc || rino: rmresf , reml		/*this is the final model that generates last panel of tab 5*/


log on
nlcom /*
*/ (rho_het: (exp(_b[lns1_1_2:_cons])^2 + varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2)) / (exp(_b[lns1_1_2:_cons]  )^2 + varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2)+ exp(_b[lnsig_e:_cons]  )^2 )) ///
   (het_share_statio:  ((_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons]  )^2))/((exp(_b[lns1_1_2:_cons])^2 + varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2)) / (exp(_b[lns1_1_2:_cons]  )^2 + varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2)+ exp(_b[lnsig_e:_cons]  )^2 ))) ///
   (het_share_nostatio: (varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2))/((exp(_b[lns1_1_2:_cons])^2 + varX*(_b[rmres:rmresf]^2+ exp(_b[lns1_1_1:_cons])^2))) )

log off
   
*****************************
*simple model under normality   
******************************   


clear all
set mem 5000m
set more off			



forvalues i=0(1)2 {
	use earnings_long,clear
	quietly summ cohort if bob==`i'
	global min_`i'=_result(5)
	global max_`i'=_result(6)
	forvalues y=${min_`i'}(3)${max_`i'} {
		use earnings_long,clear
		keep if bob == `i'
		ta year, ge(years)
		di `y'
		ge agesq=ageC^2
		ge eta=year-yob
		sort pnr eta
		regress  logearn years* ageC agesq if cohort==`y',noc
		predict res if cohort==`y',resid
		keep pnr pnrf logearn res year yob bob  age* cohort 
		keep if cohort==`y'
		compress
		save res_`i'_`y',replace
	}
}


forvalues i=0(1)2 {
	use earnings_long,clear
	quietly summ cohort if bob==`i'
	global min_`i'=_result(5)
	global max_`i'=_result(6)
	use res_`i'_${min_`i'}
	local j = ${min_`i'} + 3
	forvalues y = `j'(3)${max_`i'} {			/*puts residualised earnings back together by bob*/
		append using res_`i'_`y'
		erase res_`i'_`y'.dta
	}
	erase res_`i'_${min_`i'}.dta

	bysort pnr: ge mres=sum(res)/_N						/*long term earnings as the average earnings in that range*/
	bysort pnr: keep if _n==_N
	
	reg mres
	predict rmres`i', resid						/*take deviations from generational means*/
	keep pnrf rmres`i'    
	save res_`i',replace
}


use res_0
merge 1:1 pnrf using res_1
drop _m
merge 1:1 pnrf using res_2
drop _m

qui reg rmres0 rmres1 //ancillary regression 

   
   
 cmp ( rmres0 = ) ( rmres1 = ) ( rmres2 = )  , ind(1 1 1) 

log on
nlcom (c13: tanh(_b[atanhrho_13:_cons])*sqrt(exp(_b[lnsig_3 :_cons])*exp(_b[lnsig_1 :_cons]))) ///
	  (c12: tanh(_b[atanhrho_12:_cons])*sqrt(exp(_b[lnsig_2 :_cons])*exp(_b[lnsig_1 :_cons]))) ///
	  (c23: tanh(_b[atanhrho_23:_cons])*sqrt(exp(_b[lnsig_3 :_cons])*exp(_b[lnsig_2 :_cons]))) , post

	  
nlcom (shareIG_MVN	: (.5*_b[c12]+.5*_b[c13])/ _b[c23]) 
log off
log close

   
